****
****	county_script.do
****
****
****	Replication file for
****
****	"Health and Voting in Rural America"
****
****	Katie Cahill and Christopher Ojeda
****
****	Frontiers in Political Science (2021)
****

use "county_data.dta", clear

**
** Table 2 
sum turnout physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercent householdpoverty socialcapital totrate _differenceabsolute code_2013 covi6 if code_2013 != .


**
** Table 5

* Model 1
reg turnout physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercentp householdpoverty totrate _differenceabsolutep covi6 if voterparticipation<1 & code_2013==6, robust

* Model 2
reg turnout physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercentp householdpoverty totrate _differenceabsolutep covi6 if voterparticipation<1 & code_2013<6, robust

* Model 3
reg turnout physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercentp householdpoverty totrate _differenceabsolutep code_2013 interact_ruralhealth if voterparticipation<1, robust

* Model 4
reg turnout physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercentp householdpoverty totrate _differenceabsolutep covi6 interact_healthcovi if voterparticipation<1 & code_2013==6, robust

* Model 5
reg turnout physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercentp householdpoverty totrate _differenceabsolutep covi6 interact_healthcovi if voterparticipation<1 & code_2013<6, robust

* Model 6
reg turnout physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercentp householdpoverty totrate _differenceabsolutep code_2013 covi6 interact_ruralhealth if voterparticipation<1, robust



**
** Figure 1

**
** Figure 1
regress voterparticipation_percent c.physunhealthydays_permonth##c.code_2013 median_age highschoolorless_percent cvap_est_whitepercent householdpoverty totrate _differenceabsolute covi6 if voterparticipation_percent<1, robust
margins, at(physunhealthydays_permonth=(2.0(.25)6.5) code_2013=(1)) atmean
matrix B = r(table)'
svmat B, name(model3_code1)

regress voterparticipation_percent c.physunhealthydays_permonth##c.code_2013 median_age highschoolorless_percent cvap_est_whitepercent householdpoverty totrate _differenceabsolute covi6 if voterparticipation_percent<1, robust
margins, at(physunhealthydays_permonth=(2.0(.25)6.5) code_2013=(2)) atmean
matrix B = r(table)'
svmat B, name(model3_code2)

regress voterparticipation_percent c.physunhealthydays_permonth##c.code_2013 median_age highschoolorless_percent cvap_est_whitepercent householdpoverty totrate _differenceabsolute covi6 if voterparticipation_percent<1, robust
margins, at(physunhealthydays_permonth=(2.0(.25)6.5) code_2013=(3)) atmean
matrix B = r(table)'
svmat B, name(model3_code3)

regress voterparticipation_percent c.physunhealthydays_permonth##c.code_2013 median_age highschoolorless_percent cvap_est_whitepercent householdpoverty totrate _differenceabsolute covi6 if voterparticipation_percent<1, robust
margins, at(physunhealthydays_permonth=(2.0(.25)6.5) code_2013=(4)) atmean
matrix B = r(table)'
svmat B, name(model3_code4)

regress voterparticipation_percent c.physunhealthydays_permonth##c.code_2013 median_age highschoolorless_percent cvap_est_whitepercent householdpoverty totrate _differenceabsolute covi6 if voterparticipation_percent<1, robust
margins, at(physunhealthydays_permonth=(2.0(.25)6.5) code_2013=(5)) atmean
matrix B = r(table)'
svmat B, name(model3_code5)

regress voterparticipation_percent c.physunhealthydays_permonth##c.code_2013 median_age highschoolorless_percent cvap_est_whitepercent householdpoverty totrate _differenceabsolute covi6 if voterparticipation_percent<1, robust
margins, at(physunhealthydays_permonth=(2.0(.25)6.5) code_2013=(6)) atmean
matrix B = r(table)'
svmat B, name(model3_code6)


gen seq2a = _n
replace seq2a = . if _n > 19

forvalues i=1(1)6{
replace model3_code`i'1 = model3_code`i'1*100
replace model3_code`i'5 = model3_code`i'5*100
replace model3_code`i'6 = model3_code`i'6*100
}		
*

	
graph twoway (rspike model3_code15 model3_code16 seq2a, lcolor(black)) ///
			(scatter model3_code11 seq2a, mfcolor(black) mlcolor(black) msymbol(o)), ///
			legend(off) ///
			ytitle("Predicted Level of Turnout") ///
			ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%" 50 "50%" 60 "60%" 70 "70%" 80 "80%" 90 "90%" 100 "100%", labsize(2.5) angle(horizontal) nogrid) ///
			yline(0 10 20 30 40 50 60 70 80 90 100, lcolor(gs8)) ///
			xlabel(1 "2" 3 "2.5" 5 "3" 7 "3.5" 9 "4" 11 "4.5" 13 "5" 15 "5.5" 17 "6" 19 "6.5", labsize(2.5)) ///
			xtitle("Physical Unhealthy Days") ///
			title("Large Central Metro") ///
			yscale(titlegap(2)) ///
			xscale(titlegap(2)) ///
			scheme(s1mono) graphregion(fcolor(white)) ///
			name(model3_code1, replace)

graph twoway (rspike model3_code25 model3_code26 seq2a, lcolor(black)) ///
			(scatter model3_code21 seq2a, mfcolor(black) mlcolor(black) msymbol(o)), ///
			legend(off) ///
			ytitle("Predicted Level of Turnout") ///
			ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%" 50 "50%" 60 "60%" 70 "70%" 80 "80%" 90 "90%" 100 "100%", labsize(2.5) angle(horizontal) nogrid) ///
			yline(0 10 20 30 40 50 60 70 80 90 100, lcolor(gs8)) ///
			xlabel(1 "2" 3 "2.5" 5 "3" 7 "3.5" 9 "4" 11 "4.5" 13 "5" 15 "5.5" 17 "6" 19 "6.5", labsize(2.5)) ///
			xtitle("Physical Unhealthy Days") ///
			title("Large Fringe Metro") ///
			yscale(titlegap(2)) ///
			xscale(titlegap(2)) ///
			scheme(s1mono) graphregion(fcolor(white)) ///
			name(model3_code2, replace)

graph twoway (rspike model3_code35 model3_code36 seq2a, lcolor(black)) ///
			(scatter model3_code31 seq2a, mfcolor(black) mlcolor(black) msymbol(o)), ///
			legend(off) ///
			ytitle("Predicted Level of Turnout") ///
			ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%" 50 "50%" 60 "60%" 70 "70%" 80 "80%" 90 "90%" 100 "100%", labsize(2.5) angle(horizontal) nogrid) ///
			yline(0 10 20 30 40 50 60 70 80 90 100, lcolor(gs8)) ///
			xlabel(1 "2" 3 "2.5" 5 "3" 7 "3.5" 9 "4" 11 "4.5" 13 "5" 15 "5.5" 17 "6" 19 "6.5", labsize(2.5)) ///
			xtitle("Physical Unhealthy Days") ///
			title("Medium Metro") ///
			yscale(titlegap(2)) ///
			xscale(titlegap(2)) ///
			scheme(s1mono) graphregion(fcolor(white)) ///
			name(model3_code3, replace)		
			
graph twoway (rspike model3_code45 model3_code46 seq2a, lcolor(black)) ///
			(scatter model3_code41 seq2a, mfcolor(black) mlcolor(black) msymbol(o)), ///
			legend(off) ///
			ytitle("Predicted Level of Turnout") ///
			ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%" 50 "50%" 60 "60%" 70 "70%" 80 "80%" 90 "90%" 100 "100%", labsize(2.5) angle(horizontal) nogrid) ///
			yline(0 10 20 30 40 50 60 70 80 90 100, lcolor(gs8)) ///
			xlabel(1 "2" 3 "2.5" 5 "3" 7 "3.5" 9 "4" 11 "4.5" 13 "5" 15 "5.5" 17 "6" 19 "6.5", labsize(2.5)) ///
			xtitle("Physical Unhealthy Days") ///
			title("Small Metro") ///
			yscale(titlegap(2)) ///
			xscale(titlegap(2)) ///
			scheme(s1mono) graphregion(fcolor(white)) ///
			name(model3_code4, replace)
			
graph twoway (rspike model3_code55 model3_code56 seq2a, lcolor(black)) ///
			(scatter model3_code51 seq2a, mfcolor(black) mlcolor(black) msymbol(o)), ///
			legend(off) ///
			ytitle("Predicted Level of Turnout") ///
			ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%" 50 "50%" 60 "60%" 70 "70%" 80 "80%" 90 "90%" 100 "100%", labsize(2.5) angle(horizontal) nogrid) ///
			yline(0 10 20 30 40 50 60 70 80 90 100, lcolor(gs8)) ///
			xlabel(1 "2" 3 "2.5" 5 "3" 7 "3.5" 9 "4" 11 "4.5" 13 "5" 15 "5.5" 17 "6" 19 "6.5", labsize(2.5)) ///
			xtitle("Physical Unhealthy Days") ///
			title("Micropolitan") ///
			yscale(titlegap(2)) ///
			xscale(titlegap(2)) ///
			scheme(s1mono) graphregion(fcolor(white)) ///
			name(model3_code5, replace)			

graph twoway (rspike model3_code65 model3_code66 seq2a, lcolor(black)) ///
			(scatter model3_code61 seq2a, mfcolor(black) mlcolor(black) msymbol(o)), ///
			legend(off) ///
			ytitle("Predicted Level of Turnout") ///
			ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%" 50 "50%" 60 "60%" 70 "70%" 80 "80%" 90 "90%" 100 "100%", labsize(2.5) angle(horizontal) nogrid) ///
			yline(0 10 20 30 40 50 60 70 80 90 100, lcolor(gs8)) ///
			xlabel(1 "2" 3 "2.5" 5 "3" 7 "3.5" 9 "4" 11 "4.5" 13 "5" 15 "5.5" 17 "6" 19 "6.5", labsize(2.5)) ///
			xtitle("Physical Unhealthy Days") ///
			title("Non-Core") ///
			yscale(titlegap(2)) ///
			xscale(titlegap(2)) ///
			scheme(s1mono) graphregion(fcolor(white)) ///
			name(model3_code6, replace)				

graph combine model3_code1 model3_code2 model3_code3 model3_code4 model3_code5 model3_code6, ycommon rows(2) graphregion(fcolor(white))			
			



**
** Table 6

**
** Predicted probabilities
sum physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercent householdpoverty totrate _differenceabsolute covi6 if code_2013<5, detail
sum physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercent householdpoverty totrate _differenceabsolute covi6 if code_2013>4, detail

/*
NON-RURAL
health: 2.7 / 4.8
age: 31.7 / 46.5
high school: 25.8 / 60.6
white: 43.9 / 97.2
poverty: 6.9 / 24.4
religion: 266.81 / 713.6489
competition: 0.0254 / 0.59
cost of voting: -1.3 / 1.05

RURAL
health: 2.7 / 5.3
age: 32.9 / 51
high school: 35.1 / 65.6
white: 42.6 / 97.9
poverty: 9.6 / 30.6
religion: 247.75 / 870.1389
competition: 0.0341 / 0.6886
cost of voting: -1.22 / 1.05
*/


*Non-Rural Counties
reg voterparticipation_percent physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercent householdpoverty totrate _differenceabsolute covi6 if voterparticipation_percent<1 & code_2013<5, robust
margins, at(physunhealthydays_permonth=(2.7 4.8)) atmean
margins, at(median_age=(31.7 46.5)) atmean
margins, at(highschoolorless_percent=(25.8 60.6)) atmean
margins, at(cvap_est_whitepercent=(.449 .972)) atmean
margins, at(householdpoverty=(6.9 24.4)) atmean
margins, at(totrate=(266.81 713.6489)) atmean
margins, at(_differenceabsolute=(0.0254 0.59)) atmean
margins, at(covi6=(-1.3 1.05)) atmean

* Rural Counties
reg voterparticipation_percent physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercent householdpoverty totrate _differenceabsolute covi6 if voterparticipation_percent<1 & code_2013>5, robust
margins, at(physunhealthydays_permonth=(2.7 5.3)) atmean
margins, at(median_age=(32.9 51)) atmean
margins, at(highschoolorless_percent=(35.1 65.6)) atmean
margins, at(cvap_est_whitepercent=(.426 .979)) atmean
margins, at(householdpoverty=(9.6 30.6)) atmean
margins, at(totrate=(247.75 870.1389)) atmean
margins, at(_differenceabsolute=(0.0341 0.6886)) atmean
margins, at(covi6=(-1.22 1.05)) atmean



**
** Table 7

* Model 7
reg turnout physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercentp householdpoverty if code_2013 == 6, robust

* Model 8
reg turnout physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercentp householdpoverty totrate if code_2013 == 6, robust

* Model 9
reg turnout physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercentp householdpoverty totrate _differenceabsolutep if code_2013 == 6, robust

* Model 10
reg turnout physunhealthydays_permonth median_age highschoolorless_percent cvap_est_whitepercentp householdpoverty totrate _differenceabsolutep covi6 if code_2013 == 6, robust



**
** Figure 2

sum covi6, detail

reg voterparticipation_percent c.physunhealthydays_permonth##c.covi6 median_age highschoolorless_percent cvap_est_whitepercent householdpoverty totrate _differenceabsolute if voterparticipation_percent<1 & code_2013<5, robust
margins, at(physunhealthydays_permonth=(2.0(.1)6.5) covi6=(-1.3 .42 1.05)) atmean
matrix B = r(table)'
svmat B, name(nonrural_covi)

reg voterparticipation_percent c.physunhealthydays_permonth##c.covi6 median_age highschoolorless_percent cvap_est_whitepercent householdpoverty totrate _differenceabsolute if voterparticipation_percent<1 & code_2013>5, robust
margins, at(physunhealthydays_permonth=(2.0(.1)6.5) covi6=(-1.3 .42 1.05)) atmean
matrix B = r(table)'
svmat B, name(rural_covi)

egen seq1 = fill(1 2 3 1 2 3)
replace seq1 = . if _n > 138

egen seq2 = seq(), f(1) t(46) b(3)
replace seq2 = . if _n > 138

replace rural_covi1 = rural_covi1*100
replace rural_covi5 = rural_covi5*100
replace rural_covi6 = rural_covi6*100
replace nonrural_covi1 = nonrural_covi1*100
replace nonrural_covi5 = nonrural_covi5*100
replace nonrural_covi6 = nonrural_covi6*100

graph twoway (rspike rural_covi5 rural_covi6 seq2 if seq1 == 1, lcolor(black)) ///
			(rspike rural_covi5 rural_covi6 seq2 if seq1 == 3, lcolor(black)) ///
			(scatter rural_covi1 seq2 if seq1 == 1, mfcolor(black) mlcolor(black) msymbol(o)) ///
			(scatter rural_covi1 seq2 if seq1 == 3, mfcolor(white) mlcolor(black) msymbol(o)), ///
			legend(order(3 4) label(3 "Low Cost of Voting") label(4 "High Cost of Voting")) ///
			ytitle("Predicted Level of Turnout") ///
			ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%" 50 "50%" 60 "60%" 70 "70%" 80 "80%" 90 "90%" 100 "100%", labsize(2.5) angle(horizontal) nogrid) ///
			yline(0 10 20 30 40 50 60 70 80 90 100, lcolor(gs8)) ///
			xlabel(1 "2" 6 "2.5" 11 "3" 16 "3.5" 21 "4" 26 "4.5" 31 "5" 36 "5.5" 41 "6" 45 "6.5", labsize(2.5)) ///
			xtitle("Physical Unhealthy Days") ///
			title("Rural Counties") ///
			yscale(titlegap(3)) ///
			xscale(titlegap(2)) ///
			scheme(s1mono) graphregion(fcolor(white)) ///
			name(rural_covi, replace)

graph twoway (rspike nonrural_covi5 nonrural_covi6 seq2 if seq1 == 1, lcolor(black)) ///
			(rspike nonrural_covi5 nonrural_covi6 seq2 if seq1 == 3, lcolor(black)) ///
			(scatter nonrural_covi1 seq2 if seq1 == 1, mfcolor(black) mlcolor(black) msymbol(o)) ///
			(scatter nonrural_covi1 seq2 if seq1 == 3, mfcolor(white) mlcolor(black) msymbol(o)), ///
			legend(order(3 4) label(3 "Low Cost of Voting") label(4 "High Cost of Voting")) ///
			ytitle("Predicted Level of Turnout") ///
			ylabel(0 "0%" 10 "10%" 20 "20%" 30 "30%" 40 "40%" 50 "50%" 60 "60%" 70 "70%" 80 "80%" 90 "90%" 100 "100%", labsize(2.5) angle(horizontal) nogrid) ///
			yline(0 10 20 30 40 50 60 70 80 90 100, lcolor(gs8)) ///
			xlabel(1 "2" 6 "2.5" 11 "3" 16 "3.5" 21 "4" 26 "4.5" 31 "5" 36 "5.5" 41 "6" 45 "6.5", labsize(2.5)) ///
			xtitle("Physical Unhealthy Days") ///
			title("Non-Rural Counties") ///
			yscale(titlegap(3)) ///
			xscale(titlegap(2)) ///
			scheme(s1mono) graphregion(fcolor(white)) ///
			name(nonrural_covi, replace)

graph combine nonrural_covi rural_covi, ycommon row(1) graphregion(fcolor(white))
